#ifdef _TEST_D_
#include "Spline.h"

double C_f(double x)
{
	return 1 / (1 + x * x);
}

// f' = -2x / (1+x^2)^2

void B3_spline(int N)
{
	INTERPC interpC;
	for (int i = 0; i < N; i++)
	{
		double x = -5 + i;
		interpC.points.push_back(x);
		interpC.interpSites[x] = {C_f(x)};
	}

	interpC.interpSites[-5].push_back(10.0 / 26 / 26);
	interpC.interpSites[5].push_back(-10.0 / 26 / 26);

	Spline<1, 3, ppForm> spline = interpolate<3, cardinalB>(interpC, complete);
	vector<double> x = {-3.5, -3, -0.5, 0, 0.5, 3, 3.5};
	cout << "error: ";
	for (int i = 0; i < x.size(); i++)
	{
		cout << fabs(C_f(x[i]) - spline.evaluate(x[i])[0]) << ", ";
	}
	cout << endl;
}

void B2_spline(int N)
{
	INTERPC interpC;
	interpC.points.push_back(-5);
	interpC.interpSites[-5] = {C_f(-5)};

	for (int i = 0; i < N; i++)
	{
		double x = -9.0 / 2 + i;
		interpC.points.push_back(x);
		interpC.interpSites[x] = {C_f(x)};
	}

	interpC.points.push_back(5);
	interpC.interpSites[5] = {C_f(5)};

	Spline<1, 2, ppForm> spline = interpolate<2, cardinalB>(interpC, notAknot);

	vector<double> x = {-3.5, -3, -0.5, 0, 0.5, 3, 3.5};
	cout << "error: ";
	for (int i = 0; i < x.size(); i++)
	{
		cout << fabs(C_f(x[i]) - spline.evaluate(x[i])[0]) << ", ";
	}
	cout << endl;
}

int main()
{
	B3_spline(11);
	B2_spline(10);

	return 0;
}

#endif
